/**
 * @file src/Zuker.cpp
 *
 * @date 28-aug-2013
 *
 * @author Youri Hoogstrate
 * @author Lisa Yu
 *
 * @section LICENSE
 * yh-kt-fold can predict RNA 2D structures including K-turns.
 * Copyright (C) 2012-2013 Youri Hoogstrate
 *
 * This file is part of yh-kt-fold.
 *
 * yh-kt-fold is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * yh-kt-fold is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program. If not, see <http://www.gnu.org/licenses/>.
 */



#include "TracebackMotif.h"
#include "DotBracket.h"
#include "KTurnMotif.h"
#include "ReadData.h"
#include "Zuker.h"



/* All the following stuff should be moved into Zuker & DotBracket class */
int stack_a[MAXLEN][2];
int stacktop = -1;


short pathmatrix[MAXLEN][MAXLEN];
bool  pathmatrix_corrected_from[MAXLEN][MAXLEN];						// Set all to false
short pathmatrix_corrected_values[MAXLEN][MAXLEN];						// Set all as 0
short loopmatrix[MAXLEN][MAXLEN][2];
char  typematrix[MAXLEN][MAXLEN];										// Development stuff, verbose stuff?

static int maxpen = 3.0;												// second number in miscloop.dat -> should be moved to ReadData.cpp



/**
 * @brief Constructs /initializes the Zuker class: include parameters and init an empty dotbracket output.
 *
 * @date 05-nov-2012
 *
 * @todo move this to this->init(); and run this->init(); or rename it to this->reset();
 */
//Zuker::Zuker(parameter_hairpinlen,parameter_enableMotifFunctionality,parameter_sequence) << USE REFERENCES INSTEAD OF POINTERS TO MAKE SURE THESE ARGUMENTS ARE NEVER EMPTY !!
Zuker::Zuker()
{
	data       = new ReadData;											// Rename this to thermodynamic_data once;
	dotBracket = new DotBracket;
	
	//set_sequence(parameter_sequence)
	//set_hairpinlen(parameter_Hairpinlen)
	//set_enable_motif_functionality()
}



/**
 * @brief Destructs the class by removing nested objects.
 *
 * @date 14-may-2013
 */
Zuker::~Zuker()
{
	this->unsetSequence();
	
	delete this->data;
}



/**
 * @brief Initiator of the Zuker class; it should reset the class entirely such that all parameters are set for a new prediction.
 *
 * @date 10-may-2013
 */
void Zuker::init()
{
	this->rnaSequenceLength = 0;
	this->rnaSequence[0] = '\0';
	
	this->energy_init();
}



/**
 * @brief (re-)inits energy and traceback settings
 *
 * @date 10-may-2013
 *
 * @todo Check whether initialization has to be done until MAXLEN or until this->rnasequenceLength; probably the latter!
 */
void Zuker::energy_init()
{
	int i, j;
	//pathmatrix_corrected_from[-1][-1] = false;
	
	for(i = -1; i < MAXLEN; i++)
	{
		for(j = -1; j < MAXLEN; j++)									// j < i ??? - or preferably -> j < i-1????
		{
			pathmatrix[i][j] = NOT_YET_CALCULATED;
			pathmatrix_corrected_from[i][j] = false;
			pathmatrix_corrected_values[i][j] = NOT_YET_CALCULATED;
			
			this->Mij[i][j] = NULL;
		}
	}
}



/**
 * @brief Set the sequence in the Zuker object & reset parameters.
 *
 * @date 14-may-2013
 *
 * @param seq RNA Sequence
 *
 * @todo rename init(); to reset();
 */
void Zuker::setSequence(char arg_sequence[])
{
	init();// Rename to reset() or resetSequence() or resetZuker;
	
	int i, j;
	char tmp_nucleotide;
	
	j = 0;
	
	this->rnaSequenceLength = (unsigned) strlen(arg_sequence);
	this->rnaSequenceEnumerated = (int *) malloc(sizeof(int) * (strlen(arg_sequence) + 1));/** @todo change unsigned CHAR one day */
	
	for(i = 0; i < this->rnaSequenceLength; i++)
	{
		tmp_nucleotide = enumerateNucleotide(arg_sequence[i]);
		this->rnaSequenceEnumerated[j] = (int) tmp_nucleotide;
		
		switch(tmp_nucleotide)											// Cleans-up dirty chars like spaces, or other weird stuff
		{
			case NUCLEOTIDE_A:
				this->rnaSequence[j++] = 'A';
				break;
			case NUCLEOTIDE_C:
				this->rnaSequence[j++] = 'C';
				break;
			case NUCLEOTIDE_U:
				this->rnaSequence[j++] = 'U';
				break;
			case NUCLEOTIDE_G:
				this->rnaSequence[j++] = 'G';
				break;
			default:
				fprintf(stderr, "\nERROR: Unknown nucleotide: [%c] (char %i).\n", arg_sequence[i], arg_sequence[i]);
				exit(EXIT_FAILURE);
				break;
		}
	}
	
	this->rnaSequence[j] = '\0';
	this->rnaSequenceEnumerated[j] = NUCLEOTIDE_NONE;
	this->rnaSequenceLength = unsigned(strlen(this->rnaSequence));
	
	this->rna_sequence = this->rnaSequence;
}



/**
 * @brief Removes all variables that have been generated for a Zuker prediction.
 *
 * @date 14-may-2013
 */
void Zuker::unsetSequence()
{
	delete this->dotBracket;
	
	if(this->rnaSequenceEnumerated != NULL)
	{
		free(this->rnaSequenceEnumerated);
	}
}



/**
 * @brief Sets the minimum required length for a hairpin loop.
 *
 * @date 09-apr-2013
 *
 * @param arg_hairpinlen Mimumum length for a hairpin; default 3.
 */
void Zuker::setHairpinlen(int arg_hairpinlen = 3)
{
	this->MINIMAL_HAIRPINLOOP_SIZE = arg_hairpinlen;
}



/**
 * @brief Finds the Gibbs free energy for a stacking structure element.
 *
 * 5': i i+1 ...
 * 3': j j-1 ...
 *
 * @date 14-apr-2013
 *
 * @param arg_i Nucleotide position of the sequence, paired to j, where i < j
 * @param arg_j Nucleotide position of the sequence, paired to i, where i < j
 *
 * @todo replace '3' by min(this->MINIMAL_HAIRPINLOOP_SIZE,0)?
*/
float Zuker::energyStack(int arg_i, int arg_j)
{
	if(((arg_j - arg_i - 1) < 3) || (!isCanonicalPairing(rnaSequence[arg_i], rnaSequence[arg_j], true)))
	{
		return infinity;
	}
	else
	{
		return this->data->eparam[0] + this->data->stack[this->rnaSequenceEnumerated[arg_i]][this->rnaSequenceEnumerated[arg_j]][this->rnaSequenceEnumerated[arg_i + 1]][this->rnaSequenceEnumerated[arg_j - 1]];
	}
}



/**
 * @brief 
 *
 * @date 04-sept-2013
*/
float Zuker::energy_motif_surrounding_stack(char arg_nucleotide5p, char arg_nucleotide3p)
{
	float energy = infinity;
	
	switch(enumeratePairing(arg_nucleotide5p, arg_nucleotide3p))
	{
		case PAIRING_GC:
			energy = -2.55;
			break;
		case PAIRING_CG:
			energy = -2.10;
			break;
			
		case PAIRING_UA:
			energy += -1.50;
			break;
		case PAIRING_AU:
			energy = -1.38;
			break;
			
		case PAIRING_GU:
			energy = -1.01;
			break;
		case PAIRING_UG:
			energy += -0.78;
			break;
	}
	
	return energy;
}







/**
 * @brief Calculates a subsequence from the RNA sequence.
 *
 * Find a subsequence of the RNA sequence at position 'start' from length
 * 'len' in either the forward (forward=true) or reverse direction and
 * put it in target[].
 *
 * @date 14-may-2012
 *
 * @param start
 * @param len
 * @param target
 * @param forwards
 */
void Zuker::subSequence(int start, int len, char target[], bool forwards)
{
	int i;
	
	if(forwards)														// start -> start+1 .. start+i .. -> start+len
	{
		for(i = 0; i < len; i++)
		{
			target[i] = this->rnaSequence[i + start];
		}
	}
	else																// start+len -> start+len-1.. start+len-i .. -> start+len-len
	{
		for(int i = 0; i < len; i++)
		{
			target[i] = this->rnaSequence[start + len - i - 1];
		}
	}
	
	target[len] = '\0';
}



/**
 * @brief Finds the minimum amount of Gibbs free energy for K-Turns and provides the 2D structure in the arg_TracebackMotif argument.
 *
 * @section DESCRIPTION
 * Finds the mimimum amount of Gibbs free energy if a K-turn is found
 * -> find Gibbs energy for stacking pairing (i,j) AND for the K-turn
 *
 * @date 13-may-2013
 *
 * @param arg_i (i) Nucleotide position of the sequence, paired to j, where i < i' (ip) < j' (jp) < j
 * @param arg_j (j) Nucleotide position of the sequence, paired to i, where i < i' (ip) < j' (jp) < j
 * @param arg_ip (i') Nucleotide position of the sequence, paired to j', where i < i' (ip) < j' (jp) < j
 * @param arg_jp (j') Nucleotide position of the sequence, paired to i', where i < i' (ip) < j' (jp) < j
 *
 * @param arg_TracebackMotif If a K-Turn is found, the traceback route for its 2D structure will be put in a this variable; NULL otherwise.
 *
 * @return The minimum amount of Gibbs free energy for a K-turn element.
 *
 * @todo improve forward backward, check for both possibilities? -> backwards and forwards motif have to be in the list BOTH
*/
float Zuker::energyKTurn(int arg_i, int arg_j, int arg_ip, int arg_jp, TracebackMotif *arg_TracebackMotif, TracebackMotif **arg_TracebackMotif_p)
{
	int tmp_lenA, tmp_lenB, tmp_lenL, tmp_len5p, tmp_len3p;
	
	arg_TracebackMotif_p = (TracebackMotif **) malloc(sizeof(TracebackMotif *));
	
	float energy, tmp_energy;
	int i, j;
	
	bool forwards;
	
	tmp_lenA = (arg_ip - arg_i) - 1;									// Number of MISMATCHES between i...i'
	tmp_lenB = (arg_j - arg_jp) - 1;									// Number of MISMATCHES between j'...j
	tmp_lenL = (arg_jp - arg_ip) - 1;									// length of the loop;
	
	if(tmp_lenA < 2 || tmp_lenA > 16 || tmp_lenB < 2 || tmp_lenB > 16 || tmp_lenL < this->MINIMAL_HAIRPINLOOP_SIZE)
	{
		return infinity;
	}
	
	
	// Define forwards / backwards
	if(tmp_lenA < tmp_lenB)												// Backwards
	{
		forwards = false;
		
		tmp_len5p = (arg_j - arg_jp) - 1;
		tmp_len3p = (arg_ip - arg_i) - 1;
	}
	else																// Forwards
	{
		forwards = true;
		
		tmp_len5p = (arg_ip - arg_i) - 1;
		tmp_len3p = (arg_j - arg_jp) - 1;
	}
	
	
	
	/**
	 * @todo should be solved with binary search / trie search to reduce complexity from 64 to log(number of different K-turn sizes)
	 */
	energy = infinity;
	j = -1;
	i = 0;
	
	char tmp_motif5p[16];												// Currently no K-turns of size > 20 are found or described.
	char tmp_motif3p[16];
	
	while(this->data->kturns[i] != NULL)
	{
		if(data->kturns[i]->getSize(FIVE_PRIME) == tmp_len5p && this->data->kturns[i]->getSize(THREE_PRIME) == tmp_len3p)
		{
			if(forwards)
			{
				this->subSequence(arg_i + 1 , this->data->kturns[i]->getSize(FIVE_PRIME), tmp_motif5p, forwards);
				this->subSequence(arg_jp + 1, this->data->kturns[i]->getSize(THREE_PRIME), tmp_motif3p, (!forwards));
			}
			else
			{
				this->subSequence(arg_i + 1 , this->data->kturns[i]->getSize(THREE_PRIME), tmp_motif3p, forwards);
				this->subSequence(arg_jp + 1, this->data->kturns[i]->getSize(FIVE_PRIME), tmp_motif5p, (!forwards));
			}
			
			if(data->kturns[i]->match(tmp_motif5p, tmp_motif3p))
			{
				tmp_energy = this->data->kturns[i]->getEnergy();
				if(tmp_energy < energy)
				{
					energy = tmp_energy;
					j = i;
				}
			}
		}
		i++;
	}
	
	
	
	/**
	 * @todo Move these values to ReadData.cpp
	 * @todo Express them as average functions
	 */
	// Correction for stacked pairing (i,j)
	switch(enumeratePairing(this->rnaSequence[arg_i], this->rnaSequence[arg_j]))
	{
		case PAIRING_GC:
			energy += -2.55;
			break;
		case PAIRING_CG:
			energy += -2.10;
			break;
			
		case PAIRING_UA:
			energy += -1.50;
			break;
		case PAIRING_AU:
			energy += -1.38;
			break;
			
		case PAIRING_GU:
			energy += -1.01;
			break;
		case PAIRING_UG:
			energy += -0.78;
			break;
	}
	
	if(j > -1)
	{
		//TracebackMotif *tmp_TracebackMotif;
		//tmp_TracebackMotif = new TracebackMotif();
		
		
		//printf("0: %p = %p \n", &arg_TracebackMotif, arg_TracebackMotif);
		//printf("1: %p = %p \n", &tmp_TracebackMotif, tmp_TracebackMotif);
		
		//arg_TracebackMotif = tmp_TracebackMotif;
		//arg_TracebackMotif_p = &tmp_TracebackMotif;
		
		//arg_TracebackMotif = new TracebackMotif();					/** @todo Find some way for proper deletion -> destruction of entire Mij matrix for example */
		arg_TracebackMotif->init(data->kturns[j]->motif5pBonds, this->data->kturns[j]->getSize(5), this->data->kturns[j]->getSize(3), forwards);
		//printf(" Not NULL\n");
		
		//printf("2: %p = %p \n", &arg_TracebackMotif, arg_TracebackMotif);
	}
	//else
	//{
	//	arg_TracebackMotif = NULL;
	//}
	
	
	return energy;
}



/**
 * @brief Finds a possible AU Gibbs free energy penalty for a Stacking element.
 *
 * @date 09-apr-2013
 *
 * @param arg_n1 Nucleotide position of the sequence, paired to j, where i < j.
 * @param arg_n2 Nucleotide position of the sequence, paired to i, where i < j.
 *
 * @return The amount of additional Gibbs free energy as penalty for having an AU.
 */
float Zuker::energyAuPenalty(char arg_n1, char arg_n2)
{
	int tmp_pairing = enumeratePairing(arg_n1, arg_n2);
	
	if(tmp_pairing == PAIRING_AU || tmp_pairing == PAIRING_UA)
	{
		return this->data->eparam[EPARAM_AU_PENALTY];
	}
	else
	{
		return 0;
	}
}



/**
 * @brief Finds the Gibbs free energy for a hairpin loop structure element.
 *
 * @section DESCRIPTION
 * The hairpinloops Gibbs free energy levels are from a fixed table upto
 * a length of 30 unpaired nucleotides. Extrapoliation is applied for a
 * number of unpaired nucleotides which >= 30.
 *
 * 5': i i+1 ...
 * 3': j j-1 ...
 *
 * @date 14-may-2012
 *
 * @param i Nucleotide position of the sequence, paired to j, where i < j.
 * @param j Nucleotide position of the sequence, paired to i, where i < j.
 *
 * @return Amount of Gibbs free energy if (i,j) enclose a hairpin loop.
 *
 * @todo Move Gibbs free energy values to 'src/ReadData.cpp'.
 */
float Zuker::energyHairpinLoop(int i, int j)
{
	int nUnpaired, key, k;
	float energy;
	nUnpaired = j - i - 1;												// The number of unpaired nucleotides between i and j: if i = j+1, unPaired = 0
	
	if(nUnpaired < this->MINIMAL_HAIRPINLOOP_SIZE)						// Too small loop
	{
		energy = infinity;
	}
	else
	{
		if(nUnpaired <= 30)
		{
			energy = this->data->hairpin[nUnpaired];
		}
		else
		{
			// Delta G = Delta G[n = 30] + prelog_constant * log(n/30)
			energy = this->data->hairpin[29] + this->data->prelog * log(nUnpaired / 30);
		}
		
		
		// Terminal mismatch
		energy += this->data->tstkh[this->rnaSequenceEnumerated[i]][this->rnaSequenceEnumerated[j]][this->rnaSequenceEnumerated[i + 1]][this->rnaSequenceEnumerated[j - 1]];
		
		
		
		if(nUnpaired == 3)													// Triloop and tetraloop. - Note triloop file is empty -
		{
			energy += this->energyAuPenalty(rnaSequence[i], rnaSequence[j]);/** @todo: this->energyAuPenalty(i,j) */
		}
		else if(nUnpaired == 4)
		{
			key = (this->rnaSequenceEnumerated[i + 5] + 1) * 32768 + (this->rnaSequenceEnumerated[i + 4] + 1) * 4096 + (this->rnaSequenceEnumerated[i + 3] + 1) * 512 + (this->rnaSequenceEnumerated[i + 2] + 1) * 64 + (this->rnaSequenceEnumerated[i + 1] + 1) * 8 + (this->rnaSequenceEnumerated[i + 0] + 1);
			
			/**
			 * @todo Translate this into a while-loop
			 * @todo move MAXTLOOPS to Zuker class and Zuker.h
			 */
			for(k = 0; k < MAXTLOOPS; k++)
			{
				if(this->data->tloop[k][0] == key)
				{
					energy += this->data->tloop[key][1];
					break;
				}
			}
		}
		
		k = i + 1;
		while(rnaSequence[k] == 'C' && k < j)								// Check for a poly-C loop, remove uppercase
		{
			k++;
		}
		
		if(k == j)
		{
			if(nUnpaired == 3)
			{
				energy += this->data->eparam[13];
			}
			else
			{
				energy += this->data->eparam[12] + nUnpaired * this->data->eparam[11];
			}
		}
		
		
		
		// CCC hairpin loop
		/**
		 * @todo changethis into if seqmatch(..,..) &&  remove uppercase
		 */
		if((i >= 2) && this->rnaSequenceEnumerated[i] == NUCLEOTIDE_G && this->rnaSequenceEnumerated[i - 1] == NUCLEOTIDE_G && this->rnaSequenceEnumerated[i - 2] == NUCLEOTIDE_G && this->rnaSequenceEnumerated[j] == NUCLEOTIDE_U)
		{
			energy += this->data->eparam[EPARAM_GGG_U_PENALTY];
		}
		
		energy += this->data->eparam[EPARAM_MISCLOOP_PENALTY];									// Contribution from miscloop
	}
	
	return energy;
}



/**
 * @brief Energy function for Bulge Loop
 *
 * @section DESCRIPTION
 * This function returns the amount of energy for a bulge loop element. The
 * parameters must be: i<ip<jp<j where in particular:
 *
 * (ip-i == 1) || (j-jp == 1)
 *
 * There are two types of bulge loops:
 *
 * - Casus "5 prime bulge":
 * 5': ... i , ... , ip ...
 * 3': ... j    ,    jp ...
 *
 * - Casus "3 prime bulge":
 * 5': ... i    ,    ip ...
 * 3': ... j , ... , jp ...
 *
 * The corresponding energy is found in additional tables.
 *
 * @date 14-may-2013
 *
 * @param i Nucleotide position of the sequence, paired to j, where i < i' (ip) < j' (jp) < j
 * @param j Nucleotide position of the sequence, paired to i, where i < i' (ip) < j' (jp) < j
 * @param ip (i') Nucleotide position of the sequence, paired to j', where i < i' (ip) < j' (jp) < j
 * @param jp (j') Nucleotide position of the sequence, paired to i', where i < i' (ip) < j' (jp) < j
 *
 * @return Amount of corresponding Gibbs free energy if there is an bulge loop betweein (i,j) and (i',j'); infinity otherwise.
 */
float Zuker::energyBulgeLoop(int i, int j, int ip, int jp)
{
	int nUnpaired;
	float energy;
	
	nUnpaired = ((ip - i - 1 > j - jp - 1) ? ip - i - 1 : j - jp - 1);
	
	if(nUnpaired < 1)													// (ip-i == 1) || (j-jp == 1) Add this constraint?
	{
		energy = infinity;
	}
	else if(nUnpaired == 1)												// Size of bulge is one nucleotide
	{
		energy = this->data->bulge[nUnpaired] + this->data->eparam[1] + this->data->stack[this->rnaSequenceEnumerated[i]][this->rnaSequenceEnumerated[j]][this->rnaSequenceEnumerated[ip]][this->rnaSequenceEnumerated[jp]];
	}
	else if(nUnpaired <= 30)											// Size of bulge is less or equal than 30 nucleotides
	{
		energy = this->data->bulge[nUnpaired] + this->data->eparam[1] + energyAuPenalty(rnaSequence[i], rnaSequence[j]) + energyAuPenalty(rnaSequence[ip], rnaSequence[jp]);
	}
	else																// Size of bulge is larger than 30 nucleotides (not in dataset->estimation)
	{
		energy = this->data->bulge[29] + this->data->prelog * log(nUnpaired / 30) + this->data->eparam[1] + energyAuPenalty(rnaSequence[i], rnaSequence[j]) + energyAuPenalty(rnaSequence[ip], rnaSequence[jp]);
	}
	
	return energy;
}



/**
 * @brief Energy function for Interior Loop
 *
 * @section DESCRIPTION
 * 5': i ... ip ... : 3'
 * 3': j ... jp ...	: 5'
 *
 * Check code for:
 * - loop contains end or beginning of sequence
 * - loop internal loop < min hairpin loop?
 *
 * @date 14-may-2013
 *
 * @param i Nucleotide position of the sequence, paired to j, where i < i' (ip) < j' (jp) < j
 * @param j Nucleotide position of the sequence, paired to i, where i < i' (ip) < j' (jp) < j
 * @param ip (i') Nucleotide position of the sequence, paired to j', where i < i' (ip) < j' (jp) < j
 * @param jp (j') Nucleotide position of the sequence, paired to i', where i < i' (ip) < j' (jp) < j
 *
 * @return Amount of corresponding Gibbs free energy if there is an internal/interior loop betweein (i,j) and (i',j'); infinity otherwise.
*/
float Zuker::energyInternalLoop(int i, int j, int ip, int jp)
{
	int nUnpaired, nAsym, pairing1, pairing2;
	float energy;
	
	nUnpaired = ip - i - 1 + j - jp - 1;
	nAsym = ip - i - 1 - (j - jp - 1);
	
	energy = this->data->eparam[EPARAM_INTERIOR_LOOP_PENALTY];
	
	if(nUnpaired > 4 || abs(nAsym) > 1)
	{
		if(nUnpaired < 30)
		{
			energy += this->data->inter[nUnpaired];
		}
		else
		{
			energy += this->data->inter[29] + this->data->prelog * log(nUnpaired / 30);
			energy += this->data->tstki[this->rnaSequenceEnumerated[i]][this->rnaSequenceEnumerated[j]][this->rnaSequenceEnumerated[i + 1]][this->rnaSequenceEnumerated[j - 1]] + this->data->tstki[this->rnaSequenceEnumerated[ip]][this->rnaSequenceEnumerated[jp]][this->rnaSequenceEnumerated[ip - 1]][this->rnaSequenceEnumerated[jp + 1]] + min(maxpen, (abs(nAsym) * this->data->poppen[min(4, min(ip - i - 1, j - jp - 1))]));
		}
	}
	else
	{
		// 5': i i+1 ip
		// 3': j j-1 jp
		if (nUnpaired == 2 && nAsym == 0)
		{
			pairing1 = enumeratePairing(rnaSequence[i], rnaSequence[j]);
			pairing2 = enumeratePairing(rnaSequence[ip], rnaSequence[jp]);
			
			if(pairing1 != PAIRING_NONE && pairing2 != PAIRING_NONE)
			{
				energy += this->data->sint2[pairing1][pairing2][this->rnaSequenceEnumerated[i + 1]][this->rnaSequenceEnumerated[j - 1]];
			}
			else
			{
				energy = infinity;
			}
		}
		
		// 5': i i+1 ip
		// 3': j j-1 j-2 jp
		else if (nUnpaired == 3 && nAsym == -1 )
		{
			pairing1 = enumeratePairing(rnaSequence[i], rnaSequence[j]);
			pairing2 = enumeratePairing(rnaSequence[ip], rnaSequence[jp]);
			
			if(pairing1 != PAIRING_NONE && pairing2 != PAIRING_NONE)
			{
				energy += this->data->asint3[pairing1][pairing2][this->rnaSequenceEnumerated[i + 1]][this->rnaSequenceEnumerated[j - 1]][this->rnaSequenceEnumerated[j - 2]];
			}
			else
			{
				energy = infinity;
			}
		}
		
		// 5': i i+1 i+2 ip jp jp+1 j :3'
		// rotate 180 degree
		// 3': j j-1 jp ip ip-1 ip-2 i :5'
		else if(nUnpaired == 3 && nAsym == 1)
		{
			pairing1 = enumeratePairing(rnaSequence[jp], rnaSequence[ip]);
			pairing2 = enumeratePairing(rnaSequence[j], rnaSequence[i]);
			
			if(pairing1 != PAIRING_NONE && pairing2 != PAIRING_NONE)
			{
				energy += this->data->asint3[pairing1][pairing2][this->rnaSequenceEnumerated[jp + 1]][this->rnaSequenceEnumerated[ip - 1]][this->rnaSequenceEnumerated[ip - 2]];
			}
			else
			{
				energy = infinity;
			}
		}
		
		// 5': i i+1 i+2 ip :3'
		// 3': j j-1 j-2 jp :5'
		else if(nUnpaired == 4 && nAsym == 0)
		{
			pairing1 = enumeratePairing(rnaSequence[i], rnaSequence[j]);
			pairing2 = enumeratePairing(rnaSequence[ip], rnaSequence[jp]);
			
			if(pairing1 != PAIRING_NONE && pairing2 != PAIRING_NONE)
			{
				energy += this->data->sint4[pairing1][pairing2][this->rnaSequenceEnumerated[i + 1]][this->rnaSequenceEnumerated[j - 1]][this->rnaSequenceEnumerated[i + 2]][this->rnaSequenceEnumerated[j - 2]];
			}
			else
			{
				energy = infinity;
			}
		}
		else
		{
			energy = infinity;
		}
	}
	
	return energy;
}



/**
 * @brief Vij Function - energy if (i,j) pair, otherwise return infinity
 *
 * @date 10-may-2013
 *
 * @param i Nucleotide position of the sequence, paired to j
 * @param j Nucleotide position of the sequence, paired to i
 *
 * @return amount of Gibbs free energy provided for folding nucleotide i with j assuming i and j are paired
 */
float Zuker::v(int i, int j)
{
	int ip, jp;
	float energy, tmp, tmp_k;
	
	TracebackMotif *tmp_TracebackMotif = new TracebackMotif;
	TracebackMotif **tmp_TracebackMotif_p = NULL;						// Useless; for development purpose only.
	
	std::string substring_5p;
	std::string substring_3p;
	
	Motif *tmp_motif;
	
	if(pathmatrix[i][j] > NOT_YET_CALCULATED)							// already calculated; dynamic programming
	{
		return this->Vij[i][j];
	}
	
	if(!isCanonicalPairing(rnaSequence[i], rnaSequence[j], true) || (j - i) < this->MINIMAL_HAIRPINLOOP_SIZE + 1)
	{
		this->Vij[i][j] = infinity;
		return infinity;
	}
	
	energy = this->energyHairpinLoop(i, j);								// Hairpin element
	typematrix[i][j] = 'H';
	
	loopmatrix[i][j][0] = i + 1;
	loopmatrix[i][j][1] = j - 1;
	
	/** @todo Create function setMotfiObject(i,j, *MotifObject), getMotifObject(i,j), deleteMotifObject(i,j) */
	delete this->Mij[i][j];												// Reset K-turn flag
	this->Mij[i][j] = NULL;
	this->Nij[i][j] = NULL;
	
	for(ip = i + 1; ip < j; ip++)										// Chechs for every 'K' what the prefered structure is
	{
		for(jp = j - 1; jp > ip; jp--)
		{
			if(ip == (i + 1) && jp == (j - 1))							// Stacking element
			{
				tmp = this->energyStack(i, j) + this->v(i + 1, j - 1);
				if (tmp < energy)
				{
					energy = tmp;
					
					loopmatrix[i][j][0] = i + 1;
					loopmatrix[i][j][1] = j - 1;
					typematrix[i][j] = 'S';
					
					delete this->Mij[i][j];								// Reset K-turn flag
					this->Mij[i][j] = NULL;
					
					
					this->Nij[i][j] = NULL;
				}
			}
			else if(ip == (i + 1) || jp == (j - 1))						// Bulge-loop element
			{
				tmp = this->energyBulgeLoop(i, j, ip, jp) + this->v(ip, jp);
				if(tmp < energy)
				{
					energy = tmp;
					
					loopmatrix[i][j][0] = ip;
					loopmatrix[i][j][1] = jp;
					typematrix[i][j] = 'B';
					
					delete this->Mij[i][j];								// Reset K-turn flag
					this->Mij[i][j] = NULL;
					
					
					this->Nij[i][j] = NULL;
				}
			}
			else														// Interior loop or K-turn element
			{
				tmp = this->energyInternalLoop(i, j, ip, jp) + this->v(ip, jp);
				
				// Kills a lot of time: !!!
				tmp_k = this->energyKTurn(i, j, ip, jp, tmp_TracebackMotif, tmp_TracebackMotif_p) + this->v(ip, jp);
				
				/** @note FUTURISTIC */
				std::string substring_5p = this->rna_sequence.substr(i+1,ip-i-1);
				std::string substring_3p = this->rna_sequence.substr(jp+1,j-jp-1);
				
				tmp_motif = this->data->motifs.search_motif( substring_5p ,substring_3p );
				
				if(tmp_motif != NULL)
				{
					//tmp_k = tmp_motif->gibbs_free_energy + this->energy_motif_surrounding_stack(this->rna_sequence[i], this->rna_sequence[j]) + this->v(ip, jp);
				}
				else
				{
					tmp_k = infinity;
				}
				
				
				if(tmp < energy)
				{
					energy = tmp;
					
					loopmatrix[i][j][0] = ip;
					loopmatrix[i][j][1] = jp;
					typematrix[i][j] = 'I';
					
					/** @deprecated */
					delete this->Mij[i][j];								// Reset K-turn flag
					this->Mij[i][j] = NULL;
					
					/** @note FUTURISTIC */
					this->Nij[i][j] = NULL;
				}
				
				if(tmp_k < energy)										// K-turn element
				{
					/*
					./yh-kt-fold -s GGGAGUCUUGUGAUGAGAAGUACUGGAUCUGAAGUAGCCCUUUUUGGGCUACUUGUGAUGAAACACUCAUGGUCUGAAGACUCCC
					!! YES:  original: 999986.375000   new: 1999985.375000 
					[7][77] K-Turn
					 [11,76]*
					 [12,75]*
					 [13,74]*
					[28][59] K-Turn
					 [31,56]*
					 [30,57]*
					 [29,58]*
					Sequence length: 85
					 */
					 /*
					if(tmp_motif == NULL)
					{
						printf("Problem %i,%i,%i,%i [%s][%s]\n",i,ip,jp,j,substring_5p.c_str(),substring_3p.c_str());
					}
					else
					{
						printf("!! YES:  original: %f   new: %f \n",tmp_k,tmp_motif->gibbs_free_energy + this->energy_motif_surrounding_stack(this->rna_sequence[i], this->rna_sequence[j]) + this->v(ip, jp));
					}
					*/
					
					energy = tmp_k;
					
					loopmatrix[i][j][0] = ip;
					loopmatrix[i][j][1] = jp;
					typematrix[i][j] = 'K';
					
					this->Mij[i][j] = tmp_TracebackMotif;				// Set K-Turn flag.
					
					this->Nij[i][j] = tmp_motif;						// Motif element
					
					// Create new one
					tmp_TracebackMotif = new TracebackMotif();			/** @todo The creation of the traceback motif should happen in the energyKTurn() function! */
				}
			}
		}
		
		if(ip < j - 1)													// Multi-loop element
		{
			tmp = this->w(i + 1, ip) + this->w(ip + 1, j - 1);
			
			if(tmp < energy)
			{
				energy = tmp;
				
				loopmatrix[i][j][0] = ip;
				loopmatrix[i][j][1] = ip;
				typematrix[i][j] = 'V';
				
				delete this->Mij[i][j];
				this->Mij[i][j] = NULL;									// Reset K-turn flag
				
				
				this->Nij[i][j] = NULL;
			}
		}
	}
	
	delete tmp_TracebackMotif;
	
	this->Vij[i][j] = energy;
	return energy;
}



/**
 * @brief Wij Function - provided Gibbs free energy for sequence i...j
 *
 * @date 01-mar-2013
 *
 * @param i Nucleotide position of the sequence, paired to j, where i < j
 * @param j Nucleotide position of the sequence, paired to i, where i < j
 *
 * @return amount of Gibbs free energy provided for folding nucleotide i with j
 */
float Zuker::w(int i, int j)
{
	if(pathmatrix[i][j] > NOT_YET_CALCULATED)							// Already calculated; dynamic programming function
	{
		return this->Wij[i][j];
	}
	
	int k;
	float energy, tmp;
	
	typematrix[j][i] = ' ';
	
	if((j - i) < (this->MINIMAL_HAIRPINLOOP_SIZE + 1))
	{
		this->v(i, j);
		energy = 0;
		pathmatrix[i][j] = UNBOUND;
		
		pathmatrix_corrected_values[i][j] = pathmatrix[i][j];
		typematrix[i][j] = 'U';
	}
	else
	{
		energy = this->v(i, j);
		
		if(!isCanonicalPairing(rnaSequence[i], rnaSequence[j], true))
		{
			pathmatrix[i][j] = UNBOUND;
		}
		else
		{
			pathmatrix[i][j] = BOUND;
		}
		
		pathmatrix_corrected_values[i][j] = pathmatrix[i][j];
		pathmatrix_corrected_from[i][j] = true;
		
		for(k = i; k < j; k++)											// Find bifurcation in non-paired region
		{
			tmp = this->w(i, k) + this->w(k + 1, j);
			
			if(tmp < energy)
			{
				pathmatrix_corrected_from[i][j] = false;
				
				energy = tmp;
				pathmatrix[i][j] = k;									// assign loop/jump?
				typematrix[i][j] = 'W';
			}
		}
	}
	
	this->Wij[i][j] = energy;
	return energy;
}



/**
 * @brief Calculate the Gibbs free energy; fill the matrices Vij and Wij.
 *
 * @section DESCRIPTION
 * The Vij, Wij and Mij matrix are filled in this function. This is done
 * by walking over the diagonals of the matrix to avoid a recursive
 * explosion. Each position in a diagnoal can be calculated independently
 * from each other.
 *
 * @date 23-may-2013
 *
 * @todo Return: energy at i,j
 */
void Zuker::GibbsFreeEnergy()
{
	int i, j;
	
	for(j = 0; j < this->rnaSequenceLength; j++)
	{
		#pragma omp parallel for										// Paralelization / threading: "still reachable" memory error..
		for(i = 0; i < this->rnaSequenceLength - j; i++)
		{
			this->w(i, i + j);
		}
	}
}



/**
 * @brief Push (i,j) & matrix-flag onto the stack
 *
 * @date 06-dec-2012
 *
 * @param i Nucleotide position of the sequence, paired to j, where i < j
 * @param j Nucleotide position of the sequence, paired to i, where i < j
 * @param pick_from_v_path
*/
void Zuker::tracebackPush(int i, int j, bool pick_from_v_path)
{
	stacktop++;
	stack_a[stacktop][0] = i;
	stack_a[stacktop][1] = j;
	
	stack_a[stacktop][2] = pick_from_v_path;
}



/**
 * @brief Pops from the stack.
 *
 * @date 05-nov-2012
 *
 * @param i Nucleotide position of the sequence, paired to j, where i < j
 * @param j Nucleotide position of the sequence, paired to i, where i < j
 * @param pick_from_v_path
 *
 * @todo this->stacktop?
 *
 * @return True for success; False otherwise.
 */
bool Zuker::tracebackPop(int *i, int *j, bool *pick_from_v_path)
{
	if(stacktop >= 0)
	{
		(*i) = stack_a[stacktop][0];
		(*j) = stack_a[stacktop][1];
		
		(*pick_from_v_path) = stack_a[stacktop][2];
		
		stacktop--;
		
		return true;
	}
	
	return false;
}



/**
 * @brief The traceback algorithm, finds the optimal path in the matrices.
 *
 * @section DESCRIPTION
 * This function finds the path back. It will choose between
 * the route provided by V or W scoring.
 *
 * @date 15-may-2013
 */
void Zuker::traceback()
{
	int i, j, k, ip, jp;
	
	int tmp_i, tmp_j;
	
	bool pick_from_v_path;
	
	this->tracebackPush(0, this->rnaSequenceLength - 1, false);
	
	while(this->tracebackPop(&i, &j, &pick_from_v_path))
	{
		if(i < j)
		{
			k = pathmatrix[i][j];
			
			ip = loopmatrix[i][j][0];
			jp = loopmatrix[i][j][1];
			
			
			if(pick_from_v_path || pathmatrix_corrected_from[i][j])		// Decide which matrix to pick from
			{
				k = pathmatrix_corrected_values[i][j];
				pick_from_v_path = true;
			}
			
			
			if(k == BOUND)												// Parse the route
			{
				dotBracket->store(i, j);								// Eventual storage of pairings
				
				
				if(this->Nij[i][j] != NULL)								// If a Motif is found, trace internal structure also back
				{
					/*printf("[%d][%d] K-Turn [%s", i, j,this->Nij[i][j]->name.c_str());
					if(this->Nij[i][j]->is_reversed)
					{
						printf(" - REVERSED");
					}
					printf("]*\n");*/
					
					while(this->Nij[i][j]->traceback_pop(&tmp_i, &tmp_j))
					{
						//printf(" [%i,%i]\n",ip+tmp_i,jp+tmp_j);
						this->dotBracket->store(ip+tmp_i,jp+tmp_j);
					}
					
				}
				
				
				/*
				if(this->Mij[i][j] != NULL)								// If a Motif is found, trace internal structure also back
				{
					printf("[%d][%d] K-Turn\n", i, j);					// @todo Create verbose / quiet flag wether te report this or not 
					while(this->Mij[i][j]->pop(&ip, &jp))
					{
						printf(" [%i,%i]*\n",i + ip + 1, j - jp - 1);
						this->dotBracket->store(i + ip + 1, j - jp - 1);
					}
					
					delete this->Mij[i][j];								// @todo Create a separate destructor since objects that fall outsite the alignment are not deleted.
				}
				*/
				
				// Continue with enclosing base pair.
				ip = loopmatrix[i][j][0];
				jp = loopmatrix[i][j][1];
				
				if(ip != jp)
				{
					this->tracebackPush(ip, jp, pick_from_v_path);
				}
				else
				{
					this->tracebackPush(i + 1, ip, false);
					this->tracebackPush(ip + 1, j - 1, false);
				}
			}
			else if(k >= 0)
			{
				this->tracebackPush(i, k, false);
				this->tracebackPush(k + 1, j, false);
			}
		}
	}
}



/**
 * @brief Push (i,j) & matrix-flag onto the stack (as written by Lisa Yu)
 *
 * @date 15-feb-2013
 * 
 * @author Lisa Yu
 *
 * @param i Nucleotide position of the sequence, paired to j, where i < j
 * @param j Nucleotide position of the sequence, paired to i, where i < j
*/
void Zuker::tracebackOriginalPush(int i, int j)
{
	stacktop++;
	
	stack_a[stacktop][0] = i;
	stack_a[stacktop][1] = j;
}



/**
 * @brief Pops from the stack (as written by Lisa Yu)
 *
 * @date 05-feb-2013
 * 
 * @author Lisa Yu
 *
 * @param i Nucleotide position of the sequence, paired to j, where i < j
 * @param j Nucleotide position of the sequence, paired to i, where i < j
 *
 * @return success or not success
 */
bool Zuker::tracebackOriginalPop(int *i, int *j)
{
	if(stacktop >= 0)
	{
		(*i) = stack_a[stacktop][0];
		(*j) = stack_a[stacktop][1];
		
		stacktop--;
		
		return true;
	}
	
	return false;
}



/**
 * @brief The 'oroginal' traceback algorithm as written by Lisa Yu, finds the optimal path in the matrices.
 *
  * @author Lisa Yu
 * 
 * @details This function finds the path back. It will choose between
 * the route provided by V or W scoring. It's known to be sub-optimal.
 *
 * @date 15-may-2013
 */
void Zuker::tracebackOriginal()
{
	int i, j, k, ip, jp;
	
	this->tracebackOriginalPush(0, this->rnaSequenceLength - 1);
	
	while(this->tracebackOriginalPop(&i, &j))
	{
		if(i < j)
		{
			k = pathmatrix[i][j];
			
			ip = loopmatrix[i][j][0];
			jp = loopmatrix[i][j][1];
			
			
			
			// Parse route
			if(k == BOUND)
			{
				this->dotBracket->store(i, j);							// Eventual storage of pairings.
				
				ip = loopmatrix[i][j][0];								// Continue with enclosing base pair.
				jp = loopmatrix[i][j][1];
				
				if(ip != jp)
				{
					this->tracebackOriginalPush(ip, jp);
				}
				else
				{
					this->tracebackOriginalPush(i + 1, ip);
					this->tracebackOriginalPush(ip + 1, j - 1);
				}
			}
			else if(k >= 0)
			{
				this->tracebackOriginalPush(i, k);
				this->tracebackOriginalPush(k + 1, j);
			}
		}
	}
}



/**
 * @brief Prints the content of Vij (debugging function)
 *
 * @date 05-nov-2012
 *
 * @deprecated This is a DEBUGGING / DEVELOPMENT function.
 *
 * @param replaceInf
 */
void Zuker::printVij(bool replaceInf = true)
{
	int i, j;
	printf("\n\nVij:\n\n");
	
	for(i = 0; i < this->rnaSequenceLength; i++)
	{
		for(j = 0; j < this->rnaSequenceLength; j++)
		{
			if(this->Vij[i][j] != infinity || replaceInf == false)
			{
				printf("%7.2lf", this->Vij[i][j]);
			}
			else
			{
				printf("%7s\t", "inf");
			}
		}
		printf("\n");
	}
	printf("\n\n");
}



/**
 * @brief Print the content of Wij (debugging function)
 *
 * @deprecated This is a DEBUGGING / DEVELOPMENT function.
 *
 * @date 05-nov-2012
 */
void Zuker::printWij()
{
	int i, j;
	printf("\n\nWij:\n\n");
	
	for(i = 0; i < this->rnaSequenceLength; i++)
	{
		for(j = 0; j < this->rnaSequenceLength; j++)
		{
			printf("%7.2lf", this->Wij[i][j]);
		}
		printf("\n");
	}
	printf("\n\n");
}



/**
 * @brief Content of PathMatrix (debugging function)
 *
 * @deprecated This is a DEBUGGING / DEVELOPMENT function.
 *
 * @date 05-nov-2012
 */
void Zuker::printPathMatrix()
{
	/*
	 * Wij shouldn't contain inf values if I'm right...
	 * */
	int i, j;
	printf("\n\nPathmatrix:\n\n");
	
	for(i = 0; i < this->rnaSequenceLength; i++)
	{
		for(j = 0; j < this->rnaSequenceLength; j++)
		{
			printf("%5.2d", pathmatrix[i][j]);
		}
		printf("\n");
	}
	printf("\n\n");
}



/**
 * @brief Content of Types (debugging function)
 *
 * @deprecated This is a DEBUGGING / DEVELOPMENT function.
 *
 * @date 05-nov-2012
 */
void Zuker::printTypeMatrix()
{
	/*
	 * Wij shouldn't contain inf values if I'm right...
	 * */
	int i, j;
	printf("\n\nTypes:\n\tU=Unmatched\n\tH=Hairpin\n\tS=Stack\n\tW=W-type bifurc\n\tV=V-type bifurc\n\tB=Bulge\n\tI=Interior\n\n");
	
	for(i = 0; i < this->rnaSequenceLength; i++)
	{
		for(j = 0; j < this->rnaSequenceLength; j++)
		{
			printf("%c ", typematrix[i][j]);
		}
		printf("\n");
	}
	printf("\n\n");
}



/**
 * @brief Print base pairs
 *
 * @date 05-nov-2012
 *
 * @date Change this to Zuker::output(), add boolean arguments for [dotbracket][ct][energy][xml]
 */
void Zuker::print2DStructure()
{
	printf("Sequence length: %d\n", this->rnaSequenceLength);
	printf("%s\n", rnaSequence);
	
	dotBracket->format(this->rnaSequenceLength);
	
	printf("\n");
	
	printf("Zuker dE         %g kcal/mole", this->Wij[0][this->rnaSequenceLength - 1]);
	printf("\n");
}

